%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solveur primal gradient conjugé distribué de poutre 1D.
%	solveur classique avec possibilite de parallelisme du probleme de 
%	traction d'une poutre avec un Gradient conjugue  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [sol] = Solveur_primal_GCD (nombreElements,nombreDomaines, F, donnees, u0, errMax)

keq = donnees.E*donnees.S*nombreElements*nombreDomaines/donnees.L;

% creation des matrices blocs
Kii = keq*2*eye(nombreElements-1);
Kiidess = -keq*eye(nombreElements-2);
Kii(2:end,1:end-1) = Kii(2:end,1:end-1) + Kiidess;
Kii(1:end-1,2:end) = Kii(1:end-1,2:end) + Kiidess;

Kbb = keq*eye(2);

Kbi = zeros(2,nombreElements-1);
Kbi(end,end)= -keq;
Kbi(1,1)= -keq;

Vtmp = Kii\Kbi';
Ktmp = Kbi*Vtmp;

% creation du complement de Schur primal elementaire
Spelem = Kbb-Ktmp; 

% initialisation
bp=[zeros(nombreDomaines-1,1);F];

r = zeros(nombreDomaines,1);

r(1) = bp(1)+Spelem(2,2)*u0(1);

for k=1:nombreDomaines-1
    
    r(k:k+1) =r(k:k+1) + bp(k:k+1) - Spelem*u0(k:k+1);
    
end


w = r;

z = zeros(nombreDomaines,1);
z(1) = Spelem(2,2)*w(1);



for k=1:nombreDomaines-1
    
    z(k:k+1) = z(k:k+1)+ Spelem*w(k:k+1);
    
end


u = u0;
B = 0;

% resolution par gradient conjugue

for i=1:nombreDomaines+1
  a = dot(r,r)/dot(z,w); 
  u = u+a*w;
  d = dot(r,r);
  r = r - a*z;
    if( norm(r) < errMax )
        break;
    end
    
    B = dot(r,r)/d;
    w = r + B*w;

    
    z = zeros(nombreDomaines,1);
    z(1) = Spelem(2,2)*w(1);
    
    for k=1:nombreDomaines-1
        
        z(k:k+1) = z(k:k+1)+Spelem*w(k:k+1);
        
    end

end

u = [0;u]

% reconstruction de la solution
sol= zeros(nombreElements*nombreDomaines+1,1);

for i = 1:nombreDomaines
    sol((i-1)*nombreElements+1)=u(i);
    Uinter = (-Kii\Kbi')*u(i:i+1);
    
    sol((i-1)*nombreElements+2:i*nombreElements)=Uinter;
    
end
sol(end) = u(end);

plot(sol);

end
